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Abstract 



We have developed a new numerical scheme to obtain quasiequilibrium 
structures of nonaxisymmetric compact stars such as binary neutron star 
systems as well as the spacetime around those systems in general relativ- 
ity. Although, strictly speaking, there are no equilibrium states for binary 
configurations in general relativity, the timescale of orbital change due to 
gravitational wave radiation is so long compared with the orbital period that 
we can assume nonaxisymmetric systems in "quasiequilibrium" states. 



Concerning quasiequilibrium states of binary systems in general relativity, 
several investigations have been already carried out by assuming conformal 
flatness of the spatial part of the metric. However, the validity of the con- 
formally flat treatment has not been fully analyzed except for axisymmetric 
configurations. Therefore it is desirable to solve quasiequilibrium states by 
developing totally different methods from the conformally flat scheme. In this 
paper we present a new numerical scheme to solve directly the Einstein equa- 
tions for 3D configurations without assuming conformal flatness, although 
we make use of the simplified metric for the spacetime. This new formula- 
tion is the extension of the scheme which has been successfully applied for 
structures of axisymmetric rotating compact stars in general relativity. It 
is based on the integral representation of the Einstein equations by taking 
the boundary conditions at infinity into account. We have checked our nu- 
merical scheme by computing equilibrium sequences of binary polytropic star 
systems in Newtonian gravity and those of axisymmetric polytropic stars in 
general relativity. We have applied this numerical code to binary star sys- 
tems in general relativity and have succeeded in obtaining several equilibrium 
sequences of synchronously rotating binary polytropes with the polytropic 

indices N = 0.0,0.5 and 1.0. 
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I. INTRODUCTION 

One of the most important matters for theorists of relativistic astrophysics is to construct 
rehable models for quasiequihbrium configurations of binary neutron star systems. This is 
because such systems are the most promising sources of gravitational waves which can be 
observed by the gravitational wave detectors under construction, such as LIGO (USA), 
VIRGO (France/Italy), GEO (Germany/Great Britain), and TAMA (Japan). 

After the discovery of the first binary pulsar PSR1913+16 [T|], the orbital change of 
PSR1913+16 revealed that the existence of gravitational wave did successfully explain the 
observational results f^. From more detailed observations of PSR1913+16, much informa- 
tion about neutron stars such as the mass has been obtained and at the same time many 
other binary pulsar systems have been found 1^-0] . Since there are many binary systems 
composed of neutron stars/black holes, it is estimated that gravitational waves from binary 
systems will be observed directly once a year in galaxies within 200Mpc ||^. Gravitational 
wave detectors will become new eyes in the 21st century for uncovering extreme states of 
the universe. 

In order to understand coalescing stages of compact binary systems from observation of 
gravitational waves, we need to know detailed processes of coalescing or merging phases. This 
can be done by performing dynamical computations of evolution of binary star systems to a 
very high accuracy. However, it is still difficult to do highly accurate numerical computations 
of 3D configurations in general relativity. One reason for that is the very existence of 
gravitational wave radiation from the systems. There have been many attempts to take this 
approach and many important results have been obtained (see, e.g., Oohara, Nakamura and 



Shibata 0, Nakamura ||T0[); but fully satisfactory results have not been obtained yet. 

Recently several groups have begun to attack this problem from a different standpoint. 
Although binary systems emit gravitational waves, the timescale of its effect on the orbital 
motion is rather long compared with the orbital period except for the final stage of coales- 
cence because the energy loss rate or the angular momentum loss rate due to gravitational 



waves is so small except for such a final phase. The order can be estimated as O [{v/c)^] 
where v and c are the typical velocity of the system and the velocity of light, respectively. 
Thus we can assume that the systems are in definite orbits. If the binary systems are syn- 
chronously rotating, the systems can be approximately treated to be in equilibrium states 
by choosing frames which rotate with the same rotational periods as the orbital periods. In 
other words, we can treat the systems in 'quasiequilibrium'. 



From a different context, in Newtonian gravity, Hachisu and Eriguchi ITT HT^] obtained 
first numerically exact results for binary configurations of polytropes including incompress- 
ible models. In their papers, secular instability of synchronously rotating binary systems 
and dynamical instability of asynchronously rotating binary systems were discussed. On the 
other hand, recent investigations in Newtonian gravity are deeply related to the problem 
mentioned above. Lai, Rasio and Shapiro [p"5| -[T8[| have studied this problem by employing 
the ellipsoidal approximation scheme in which configurations are assumed to be exact el- 
lipsoids. More recently Uryu and Eriguchi ||T9|- p^] have developed a scheme for irrotational 
binary star systems and obtained numerically exact stationary sequences. Newly obtained 
stationary sequences have been analyzed and much information about dynamical instability 
of evolutionary sequences of polytropic binary star systems has been obtained. 

However, Newtonian configurations cannot be directly applied to realistic evolution for 
binary compact star systems. Thus some authors have studied binary configurations in 
post-Newtonian regime [p^ -|30|l- Post-Newtonian analyses show that the critical angular 
velocity for instability is increased by 10-15 % compared with that of Newtonian gravity. It 
is necessary to solve general relativistic models for binary configurations to compute exact 
critical states against instabilities. 

The first numerical results for quasiequilibrium states of binary configurations in full 



general relativity have been obtained by Wilson, Mathews and Marronetti [^]. In order 
to treat the problem tractable, they have assumed that the spatial part of the metric is 
conformally flat (the conformally flat condition (CFG), hereafter). One of their results was 



surprising: binary neutron stars become unstable and collapse into black holes individually 
prior to merging (see also Mathews and Wilson I^Sj). The same problem has been investi- 
gated by applying almost the same formulation but by using a different numerical scheme 
by Baumgarte et al. |33| (also see Cook, Shapiro and Teukolsky |3^, Baumgarte et al. [0). 
Baumgarte et al. |3^ have found that there exist quasiequilibrium sequences of binary sys- 



tems just prior to merger, i.e. that individual collapse to black holes will not occur (see also 
11,11). 



In order to check the validity of the CFC and its accuracy. Cook, Shapiro and Teukolsky 



34 have computed equilibrium configurations of axisymmetric rotating polytropes by using 



two different schemes: 1) the scheme with the CFC and 2) the KEH scheme [^. Since 



differences of physical quantities between two schemes are less than 5 %, they have con- 
cluded that the scheme with the CFC will give reasonably accurate results even for other 



situations such as binary configurations. However, its validity is not fully understood p3 
for example, it is uncertain whether results from axisymmetric configurations may be ap- 
plied to nonaxisymmetric situations or not. In this sense, it is desirable to develop different 
schemes from that with the CFC and to compare results of two or many different schemes 
for nonaxisymmetric models. 

In this paper, we present a new numerical scheme to handle quasiequilibrium states 
of nonaxisymmetric configurations in general relativity. This scheme is an extension of 
the numerical scheme for axisymmetric configurations developed by Komatsu, Eriguchi and 
Hachisu (KEH scheme) [Q . The basic idea used in |^8[ is to transform the Einstein equations 



into integral equations by using the Green function for the Laplacian in the flat space. 
Since we can include the boundary conditions into the Green function, what we need to 
do is to solve integral equations by some means. In order to adopt the same procedure for 
nonaxisymmetric models as that for axisymmetric configurations, the crucial step is to make 
up the Laplacian in the flat space by arranging the Einstein equations. For nonaxisymmetric 
configurations, however, the Einstein equations for each metric component do not contain 
Laplacians in the fiat space even after proper arrangements. This can be done if we add 



appropriate terms to both sides of properly selected component equations of the Einstein 
equations, although the source terms of those equations contain extra derivatives in addition 
to the original sources. 

As the first step to construct realistic models for binary star systems, we have assumed 
a simplified form for the metric. The basic equations consist of combinations of dominant 
components of the Einstein equations for this simplified metric. We have successfully applied 
the above procedure to make up Laplacians in the fiat space and have succeeded in developing 
a new reliable numerical scheme. 

This paper is organized as follows. In section II, we introduce assumptions and basic 
formulation of the problem. The choice of the simplified metric and the Einstein equations 
are explained. In section III, we describe the numerical solving method briefiy. In section 

IV, results of our numerical computations are presented. We have tested our new method by 
comparing our results with (1) those of the axisymmetric configurations of fully relativistic 
and rapidly rotating polytropes and (2) those of the Newtonian binary systems. In section 

V, we discuss the validity of the present scheme and the future prospect. 

II. ASSUMPTIONS AND BASIC EQUATIONS 

A. Assumptions 

As mentioned in Introduction, nonaxisymmetric rotating objects cannot be in equilibrium 
states in the framework of general relativity. It is because gravitational waves carry away 
the energy as well as the angular momentum from the binary star systems to infinity and so 
no conserved physical quantities can exist for them. However, for almost all stages except 
for the last few milliseconds of coalescence of binary systems, the timescale of change of 
the system due to gravitational waves, tqW) is much longer than that of the orbital period, 
Torb- Specifically, Torb/Tcw ~ {v/cY and it is less than 1 % even for v ~ 0.3c. At the same 
time, the energy radiated due to gravitational wave emission is only a few percent of the 



total energy even when two stars come very close. Therefore, gravitational waves can be 
neglected for most stages of evolution of binary star systems. 

If we neglect the effect of gravitational waves, we can choose a proper rotating frame in 
which the system is regarded as a stationary one. If we let the angular velocity of this frame 
seen from the observer at infinity be Q, the following vector in the inertial frame: 



can be regarded as a Killing vector (see e.g. Bonazzola, Frieben and Gourgoulhon ||39| , ^ 



Bonazzola, Gourgoulhon and Marck |^), where t and if are the time and the azimuthal 



coordinates, respectively. Thus we assume that nonaxisymmetric systems are in quasiequi- 
librium states. 

Assumptions we will make in this paper are as follows: 

1) We will treat a binary system which consists of two equal mass stars in a circular orbit, 

although our code can treat a single nonaxisymmetric body. 

2) The binary star system is assumed to be in a stationary state in the rotating frame with 

the angular velocity Q. In other words, we will neglect gravitational waves from the 
binary star system. 

3) Axes of spins of two stars and that of the orbital motion are parallel to each other. A 

schematic figure of the system is shown in Figure |l]. The 2;— axis is the rotational axis 
of the orbital motion. The x— axis is connecting two centers of mass of the stars. The 
intersection of these two axes is defined as the origin of the system. The y— axis is 
perpendicular to the x— and ^— axes. We will call the x — y plane as the equatorial 
plane. 

4) Spins of two stars are synchronized to the orbital motion. Each star is rigidly rotating 

with the angular velocity Q if seen from a distant place. 

5) The matter of the star is perfect fluid and a polytropic relation is assumed. 



6) The matter distribution and the spacetime are assumed to be symmetric about three 
planes: the equatorial plane, the y — z plane and the x — z plane. 

B. Metric and Einstein equations 

In addition to the assumptions mentioned above, we further assume the following form 
for the metric in the spherical coordinates [r,6,ip) (hereafter we use the units of c = G = 1): 

ds^ = -e^'^de + r^ sin^ ee^^{d^ - uodtf + e^^rfr^ + r^e^"' dO^ , (2) 

where z/, /?, u, a and a' are the metric coefficients and they are functions of r, 9 and (p. 
This choice of the metric is not the most general one for quasiequilibrium states because 
we do not take into account tr- and/or tO- components of the metric. However, since the 
purpose of this paper is to show the effectiveness of our new scheme and the tip- component is 
considered to be the most dominant one among nondiagonal components, we use the above 
form for the metric in this paper. Furthermore, for simplicity, we assume the following 
condition: 

a = a . (3) 

This form of the metric becomes exact for stationarily axisymmetric configurations. Al- 
though the assumed metric form is incomplete for quasiequilibrium 3D configurations, the 
essential and technical part of a new numerical scheme for 3D configurations can be shown 
by choosing this kind of simplified metric and will be extended to more general form of the 
metric. 

As mentioned before, the equation of state for the matter is assumed to be polytropic as 
follows: 

p = Ke^+i/^ , (4) 

where p, K, e and A^ are the pressure, a constant, the energy density and the polytropic 
index, respectively. It should be noted that this polytropic relation is slightly different from 



that used in [^. This equation of state can be rewritten as follows by using the Lane-Emden 
function A: 



p = PcX 



l+N 



e = e.A^ 



(5) 
(6) 



where pc and Ec are the maximum pressure and the maximum energy density, respectively. 
The energy-momentum tensor, T''^, for the perfect fluid is written as: 



T^^ = {e+ p)u^vP + pg^"^ , 



(7) 



where u'^ and g^^ are the four velocity of the matter and the metric, respectively. Throughout 
this paper, Greek indices run from to 3. Concerning the four velocity, since an observer 
in the rotating frame sees that the matter is static, we can obtain the following relation by 
using the condition m'^m^ = — 1: 





u^ 



VT^ 







(8) 



where v is the proper velocity of the matter: 



V = r{Q — uj) sin 



(9) 



After a lengthy but straightforward calculation of the Ricci tensor from the metric (^ 
(see Appendix), combining the obtained Einstein equations appropriately (see Appendix) 
and introducing the following two variables p and 7 by 



p = p-u 

-i = j3 + v 



(10) 

(11) 



we can obtain the following equations for the metric functions p, 7, a and uj: 



Aa = Sa{r,9,ip) , 

,. 2d 2 cot 61 (9,, n^^/o. ^ / ^ X 
r ar r^ at/ 



(12) 
(13) 
(14) 
(15) 



where 



a^ 2d 1 92 1 ^9 



92 



9^ r^ sin 6 dip"^ 



(16) 



and 



-87re2"(e+p) 
-2 



l+t;2 /1 97 1 ,97^ 

-— ^ + ^cot6'-^ 

\r or r^ do 



l-v 

g2a-2/3 /^2^ Q^ 



n (97 (9a 



r^ sin 6 \d'p'^ dip dp) dip ^ 



^ \ 167re2> 



i^ + lcot^— ' 
r 9r r^ (9^ 



2 Vdr' r^^de' , 



-2 



g2Q 2/3 /^2^ /<9z/x2 "^^^ /"^a.g dp da du d(3^ 



+ 



+ 



r^ sin 6* \dp>'^ dip dp>'^ dp> dp) dp) dp) dp) 



+ 2{n-uj?e^"''^'' i^ + (—^^ , -^^a ^ ^(9a^2 , "9^(90 dv d(3^ 






2{Vt-u)e 



+2^'^~'^(^) 
dip 



2a-2u 



dpP- dp) dip'^ dp) dp) dp) dp) dp) 
d'^oo duo da dp d'j 



+ 2-^ + 



dp)"^ dp) dp) dp) dp) 



+ 



1 



(92 



r^ sin 6 dip^ 



(pe 



7/2 ^ 



(17) 



S,{r,e,v) = e''/' 



16.e>-fi^ + lcot^|) 
\r dr r^ do J 

+ -Il67re2° - (^i^ + ^cot^— ' 
2 I V r (9r r^ (9^ 



2 Vdr^ r^^de' , 



2a 2/3 / Q2^ (^^\2 ^'^ ^ (9a ^ dp da du df3^ 
r^ sin^ 9 \ dp)^ dip dp)^ dip dip dp) dip dip 
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2{n-ujye 



2.a-2Jd'P , ,dP^ 



+ 



r, d'^a ,da,r, dp da du 83^ 



+ 2{n-to)e 



2a-2u 



^ dif'^ dip dif'^ dip' 'dip dp dtp dtp 
d'^uj duj da dp d'j 



— (2 h2-^ H -) 

dp)'^ dp dp dp dp 



+2e 



2a-2v/'^^\2 



^dp' 



g2Q 2/3 / Q2j^ ,(9z/,n d'^a .da... 



r^ sin 9 \ dp"^ dp dp"^ dp ' 



dp da du d(3 ' 

dp dp dp dp ^ 



+ 2{n-ujfe 



2^2a-2u I '^ y , (__)^ _|_ 



d'^a ,da 



r, dp da dv dd 



+ 2{n-uo)e 



dp 



2a~2u 



dp'^ dp dp'^ dp dp dp dp dp ^ 
d'^uj duo da dp ^ d'^ -^ 

1 ^~^ \ ^~ 



dp^ dp dp dp dp' 



+ 



1 



92 



r^ sin Q dp'^ 



5,(r,^,^) = -47re^"(£ + p) + 



Ida 1 ^da \dv 1 dv 

■ cot 6*^:^ + --— + — cot 9- 



r dr r^ dO r dr 

dv d(3 1 dv d(3 

dr dr r^ dO dO 



dd 



^2a-2l3 



+ 



'd'^v 



r^ sin^ 6 ydp"^ 



+ 



dv_, 
^dp' 
d'^ii 



dv_d^ 
dp dp 



,da., 
'dp' 



2a-2u 



+ 



{Vt — uj)e 

1 d^a 
r^ sin^ 6 dp"^ 



^dp"^ dp dp dp 

d'^uj duo . dp d'l , 

(2—^ H 

dp'^ dp dp dp' 



,da, 

'dp' 



— e 



2q-2!./^^n2 



'dp' 



(19) 



S^{r,9,p) = e^^^/' 



-IQne 



2a 



{e + p){n-uo) 



+ uo < — Sttc 



2a 



[1 + v'^)e + 2v^p 



\-v' 



r \ dr 2dr j 



1 n ( '^P 1 <97 ' 

V-''''^Vm-2-de, 
de' ^dd' 
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g2a-2/? /Q2 



V 



+ 



dv . 



+ {n-uYe' 
dp 



+ 3( 



dip"^ ^ dip 



[dip^ ^dip' 
dv da dv d(3 \ 



,d'^a ,da,2\ da dp c^^v dv d(3^ 
dip dip dip dip dip ^ 
,d'^a /<9a;,2x 
''d^^^d^' ' 



+ 2 , 
dip dip dip dip dip 



+ (r2 - uo)e 



2a-2u 



d'^ijj duj , da dp d^ , 



dip"^ dip dip dip dip' 



dip 



-'■^-^-■i(|)^-^<|)^. 



-4 



+ 



g2Q:-2/3 
1 ^2 , 



cue 



P+7/2N 



'9^a ,9a. 



9a ^7 



+ ( — )'- 
c^v^^ 9^9 9^9 dip 



r^ sin^ ^ 9(y9^ 

If we multiply (rsin^siny)) to both sides of Equation (p!5|), we can obtain the following 
equation: 



A(r sin 6^ sin </? ■ tue''''"''''^) = S^{r, 9, ip) 



(21) 



(20) 



where 



SJr,9,ip) = rsm6smipSu;{r,6,ip) -\ ^-^TT^i^e''^"'^'^) 

r sm u dip 



(22) 



It should be noted that in these equations we make up Laplacians in the flat space by 

1 d^ 



adding terms with derivatives 



to both sides of the equations. The rest of the 



r^ sin^ 9 dip"^ 
Einstein equations are not used in this paper. Therefore, our solutions are not numerically 

exact ones which satisfy all components of the Einstein equations. This is similar to the CFC 

approach [0,BH] . We will discuss the choice of the components of the Einstein equations in 



a later section. 



C. Boundary conditions and integral representation of the Einstein equations 



Equations (|T2|) - (|T^ and (pT]) can be regarded as Poisson equations for the corresponding 
quantities if we treat the right hand sides as source terms, although source terms contain 
unknown quantities. 
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Under our assumption of quasiequilibrium states or no emission of gravitational waves, 
we need not consider that there are gravitational waves at infinity. This implies that the 
metric at infinity must be fiat. Therefore, we can transform these equations by using the 
Green function for the Laplacian in the fiat space, l/|r — r'|, into the integral equations as 
follows: 

1 /-oo ("K f27r 1 

p = -—e~^l^ dr' dff d^'sme'r'^SJr',e',^'). ., (23) 

4:77 Jo Jo Jo \r — r I 

7 = --^e-^/2/ dr' de' d^'sme'r'^SJr',e',^')- -, (24) 

4:71 Jo Jo Jo \r — r I 

a = -— dr' dO' dip' sin 9' r'^S^ir', 9', ip')- -, (25) 

47r Jo Jo Jo |r — r I 

1 /-oo rn r2TT _ 1 

rsin^sin^cj = e^^'^+^/^M dr' d9' dip' sin 9' r'^S J r', 9', ip')- r. (26) 

471 Jo Jo Jo r — r' 



The Green function in these equations is expanded as follows: 
1 



^/„(r,r')[Pn(cos^)P„(cos^') 



n=0 



+2 Y: y^^Pri^os9)P::{cos9') cosm(^ - ^') 



(27) 



where 



1 fr'Y 

- — , for r'/r < 1 

fniry) = {\yj (28) 

— I — I , for r/r' > 1 
r' \r'J ' ' 

and Pn and P™ are the Legendre polynomials and the associated Legendre functions, re- 
spectively. The asymptotic fiatness at infinity, r -^ cxo, 

P ~ 0(l/r) , (29) 

7 ~ 0{iy) , (30) 

a ~ 0(l/r) , (31) 

to ~ 0(l/r^) , (32) 

are satisfied automatically, if the source terms behave properly. 
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D. Hydrostatic equation 

The hydrostatic equation is derived from the conservation law T''7 = as follows: 

Vp-(e+p)VlnM* = . (33) 

By using the polytropic equation (^) and the rigid rotation law, it can be integrated to the 
following equation: 

(1 + A^) In (Ke^/^ + l) + z/ + - ln(l - v^) = C , (34) 

where C is a constant of integration. 

III. METHOD OF SOLUTION 

A. Model parameters 

Once the equation of state is fixed, we need to specify two parameters to determine one 
model for a rotating equilibrium configuration: one parameter which represents the strength 
of gravity and the other for the amount of rotation. For polytropes, we have to specify one 
more parameter because the constant i^ is a free parameter. In our formulation, we choose 
the following three parameters: 1) the maximum energy density, Ec, 2) the ratio of the 
maximum pressure to the maximum energy density, k: 

K = ^ , (35) 

and 3) the ratio of the shortest distance r^i (distance from the origin to point A, see Figure 
|l]) to the largest distance r^ (distance from the origin to point B, see Figure |l]), q: 

g--. (36) 

T 

The value of the quantity k, represents the strength of gravity because /^ ~ ^, where Vg and 

R 

R are the Schwarzschild radius and the radius of the star, respectively. This is regarded as 

14 



the parameter of compactness. For post-Newtonian models, k < 0.1 and for typical neutron 
stars, K ~ 0.2 - 0.4. 

The quantity q indirectly specifies the rotation rate. It may be possible to choose alterna- 
tively the angular velocity or the angular momentum to specify the rotation rate. However, 
from the numerical computational point of view, it is much better to choose the ratio of two 



distances instead of choosing the other physical quantities |T^ . 

For polytropes, we can introduce dimensionless physical quantities by using two constants 
c and G as well as the maximum energy density Ec- Since the quantity Ec does not appear 
in the basic equations, we can treat the problem with two parameters, k and q, in addition 
to the polytropic index N. In summary, once we specify these parameters, all we have to 
do is to solve for the metric coefficients p, 7, a and u, distributions of the energy density £, 
the pressure p and the angular velocity of the system Q. 

In actual numerical computations, we further normalize the radial coordinate and the 
metric functions as follows: 



r = roJfxT , (37) 



and 



P = fxp , (38) 

7 = /x7 , (39) 

a = f^a , (40) 

UJ = ^ofx^ , (41) 



where 



C2 



- - TBi: • '^=' 



no = ^A-kGeJc^ , (43) 

and fx is a dimensionless normalization constant which is implicitly determined from 
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tb = yfxro , or r^ = 1 . (44) 

By using these normalized quantities, the hydrostatic equation ( P^ is rewritten as: 

(1 + N) In {k\ + 1) + /,!>+ 1 ln(l -v^)=C , (45) 

where 

V = V^r(f2 - f^uj) smeexp{f^0 - i>)} ■ (46) 

Here 



B. Solving scheme 

Our formulation in this paper is almost the same as that used in Komatsu, Eriguchi and 
Hachisu |^. Therefore we can adopt the HSCF method ||TB[ or the KEH method 



as 



our solving scheme. The HSCF method is the extended version of the the Self- Consistent 
Field method developed by Ostriker and Mark [^|. 



The essence of the HSCF method or the KEH method can be briefly summarized as 
follows: At the beginning of the computation, we prepare initial guesses for the metric 
potentials p, 7, a and u, the energy density e, and the angular velocity Q as well as the 
quantity f^. Substituting them into right hand sides of the integral equations (p3D-(p^), we 
obtain new values of p, 7, a and u. At this point, we must solve the angular velocity Q, a 
constant of integration C and the scale parameter f^. These quantities are obtained from 
the hydrostatic equations ( ^Sf ) at points A, B (in Figure ^ and at point of the maximum 
density. Using newly obtained p, 7, a, u, Q, C and fx, we calculate the new energy density 
e from the hydrostatic equation (^3|). 

These newly obtained values of p, 7, a, uj, e, fl and fx are used as the trial guesses in 
the next iteration cycle. We repeat this procedure until the relative differences between new 

16 



values and old ones become small enough, i.e., typically less than 10~^. At this point we 
regard these values as converged ones. 

This iteration cycle is carried out by fixing A^, k and q. After we obtain one equilibrium 
model, we solve other models by changing the axis ratio q. In this way we can obtain a 
sequence of equilibrium configurations for the same A^ and k. 

In actual computations, we have used (r x ^ x </}) = (128 x 49 x 81) grid points and we 



have taken into account the summation in the expansion of the Green function (pTl) up to 
24th term, i.e. P2^4 (^os^) cos24y9. Although the integral region for integral equations (^31) 
- (p6D extends to infinity, we have covered only a finite region (0<f<2, O<0< 7r/2, 
< (y9 < 7r/2) as was done for axisymmetric configurations in |^. In Newtonian gravity 



(k ~ 10 ^), convergence is reached after 50 or so iterations, while in relativistic models 
(k ~ 0.5), it takes 200 or more iterations until converged solutions are obtained. 

C. Physical quantities 

Since there exist no exact equilibrium or stationary states, it is impossible to define con- 
served quantities for nonaxisymmetric configurations in general relativity. However, within 
our assumption of quasiequilibrium states or no gravitational waves, we can devise to define 
the approximate gravitational mass and the approximate angular momentum of the system. 

As mentioned before, since we neglect gravitational waves at infinity, the metric at infinity 
can be approximated by the fiat spacetime. Thus we assume that the metric function behaves 
as Equations (|29|) - (|3^) . Definitions of approximate quantities can be obtained as follows 
by taking the asymptotic behavior into consideration. 

For the approximate angular momentum, from the field equation of c<j, we can derive the 
following equation by arranging terms: 



1^ fr^-r^sin^^e^'^-^^U^Ij fsin^.r^sin^^e^/^-i^") 
r'^ or \ or I r sm 6 OU \ r OU I 



467rrsin0e2-+2;3(£ + P)^ 



y2 
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+ 2e2"+'^-^ 



ijj 



' d^a da. 



2 da ( 8(3 ^ dv^ 



dip \dip dip 
B. op.... /.^........ .o .M. e,.a*„. ..„. o„ ..e ,eft ,.,. .de oUM. equation 

are converted into surface integrals and can be integrated as 



'd_ 
dr 



2 2-2 

r ■ r sm 



,dij' 
dr 



dr > sin OdOdip 



d_ 



sin 6 ■ r sin 



'^-— — (i^ > rdrdip 
r oO 



r'sin^ee'^- 



,duj 
dr 



sin OdOdip + 



sm' 



^e^'^^ 



9^ 



rdrdip 



,dto' 



^4g3/3-z. "I sin'' OdOdip 

\ / r=oo 



-6 J f sin^ OdOdip 
-IGttJ , 



2J , (9cc; 6J 
where we have used the asymptotic behavior oo ~ — - and — — ~ at r ~ oo. By using 



rp O (I'V f 



this relation, we can find 



J = / r sin ge^"+2/3 (^ + PK 2 ^^^ ^^^^^^ 
J 1 — t>^ 



1 

47r 



g2a+/3-.^ 



(9^q; da I da dj3 r,'9^\ 



r^ sin OdrdOdip . 



(48) 



Here we should note that in the asymptotic fiat region, the coefiicient of the first term of 
the expansion of the metric component uj in terms of 1/r is interpreted as the twice of the 
total angular momentum of the system. Therefore, we may define the quantity J as the 
approximate angular momentum of the nonaxisymmetric configuration. 

In the same way, we can define the approximate gravitational mass as follows: 



M 



^2a+l3+iy 



(e + p) 



1 — f ^ 



+ 2p> r'^ sin OdrdOdip 



f 2r sin Oe^''^^^ (e + p) ^^r^ sin OdrdOdip 
J 1 — v^ 



(49) 



• n , . 2J du M 

where we have used the asymptotic behavior uo ~ -^ and -^r- ~ ^7 at r ~ oo. 



r^ — ana — ^^ — 



As for the rest mass, it is natural to define as: 

Mn = / Pnu^\/^ ■ r^ sin OdrdOdip 



iO 



Pou V 



[l+p/e)^VT^^- 



-.r"^ sin OdrdOdip , 



(50) 



where on = -, t-tt? is the baryon mass density. 

IV. RESULTS 
A. Numerical tests of our new code 

We have checked our code by applying it to two cases for which solutions have been 
obtained by other methods: 1) rotating and general relativistic axisymmetric sequences of 
polytropes and 2) Newtonian binary sequences of polytropes. 

As mentioned before, our formulation becomes exact for axisymmetric configurations be- 
cause the metric form (^ is the general one for the axisymmetric space-time and there can 
be exact stationary states. Komatsu et al. |^8| computed general relativistic and axisymmet- 
ric rotating polytropes by developing a 2D code. We have computed the same equilibrium 
sequences by applying our 3D code. 

In Figures H (a) and (b), we show equilibrium sequences of models starting from rp/re = 
0.9375 (nearly spherical) to rp/r^ = 0.5625 with k = 0.0001, 0.25 and 0.4 for (a) A^ = and 
(b) N = 0.5, respectively, where r^ and r^. are the polar radius and the equatorial radius 
of axisymmetric configurations, respectively. In these figures, the squared nondimensional 
angular velocity Q^ is plotted against the nondimensional angular momentum J, where 

Three curves show our 3D computational results and discrete points denote the results of 
the 2D code. As seen from these figures, our results of 3D code agree well with those of 
Komatsu et al. ||38| to within less than 0.5%. 

As for binary sequences, we have carried out computations of Newtonian binary se- 
quences. In our 3D code, Newtonian limit is treated by choosing k, = 0.0001. We have 
computed equilibrium sequences for A^ = 0.0, 0.5, 1.0 and 1.5 from q = 0.0 (contact phase) 
to g = 0.5. Our results are compared with those of Hachisu |jl3| in Figure ^ 
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As seen from this figure, except for A^ = 0.0 models, corresponding models of two different 
codes agree well each other to within less than 0.5%. For the A^ = 0.0 sequence, some of 
our models are not in good agreement with those of Newtonian computations. This may be 
because our 3D code does not treat the surface region of the star carefully for the constant 
density distributions of A^ = polytropes due to small mesh numbers within the star. In 
our computations, the star is covered with 64 grid points at most in the r— direction and its 
number may not be enough to treat drastically changing density distributions. Although it 
is possible to increase grid points enough to manage a stiff distribution of the density, we 
did not try it because such a stiff equation of state is not suitable for real neutron stars. 

B. Synchronously rotating relativistic binary sequences 

In real evolution of compact binary systems, the matter cannot be represented by a simple 
polytrope and structures of compact stars will change during evolution so that the value of 
K must be changing. Such realistic models can be treated by using the realistic equation of 
state and following quasiequilibrium sequences with constant baryon mass models. However, 
in this paper, we report only the results of binary sequences with fixed values of A^ and 
K because the main purpose of this paper is to show a new numerical scheme to treat 
nonaxisymmetric and general relativistic configurations in quasiequilibrium states. 

In Figures ^ (a), (b) and (c), we show the results of synchronously rotating binary 
equilibrium sequences from q = 0.0 to g = 0.5 with k = 0.0001, 0.1, 0.3 and 0.5, for 
(a) A^ = 0.0, (b) A^ = 0.5 and (c) A^ = 1.0 polytropes, respectively, (see also Tables in 
Appendix). 

From these figures, we can see that as the strength of gravity is increased, sequences 
are shifted to the part with a larger value of the dimensionless angular momentum, i.e. 
at the righter parts of the panels. Similar tendency can be found for Newtonian binary 
sequences such as Figure |^. In the Newtonian models, models with larger values of A^ are 
shifted to the righter part of the panel because of the density concentration to the central 
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part of the star. As discussed in [^, the effect of general relativity is to increase the mass 
concentration to the central part of the star. In other words, models with large k can be 
regarded as models with higher "effective" polytropic indices. Consequently sequences of 
larger k models are located at the part with a larger value of the dimensionless angular 
momentum in the figure. This situation can be clearly seen in Figure ^ In this figure, equi- 
density contours are shown for N = 0.5 polytropes with k, = 0.0001 (Newtonian model) 
and K = 0.51 (relativistic model). For the same value of A^, the density concentrates to the 
central region of the stars for relativistic models more than that for Newtonian models. 

The other characteristic feature of these figures is that for A^ = sequences there are 
turning points where the value of the dimensionless angular momentum becomes minimum 
but that for A^ = 0.5 and A^ = 1.0 sequences there are no turning points except for Newtonian 
case with A^ = 0.5. 

As discussed in Introduction, in the framework of Newtonian gravity, many sequences 
have been investigated. In particular, binary sequences are known to be connected to ax- 



isymmetric ones by way of Jacobi-Dumb-bell shaped sequences |^|. We show the whole 



relation of several equilibrium sequences of A^ = polytropes, i.e. from Newtonian to rela- 
tivistic and from axisymmetric to binary sequences as well as the Jacobi-Dumb-bell shaped 
sequence for Newtonian gravity in Figure §. 

V. DISCUSSION AND CONCLUSIONS 

A. Discussion 

As discussed before, the metric form (|]) is too simplified and may not be suitable for 
3D quasiequilibrium configurations. In general, the number of independent metric functions 
for stationary spacetime can be reduced to six. Thus, the most crucial point may be that 
we have not included nondiagonal components except for the t(y9-component. The results 
obtained in this paper may be affected by taking into account such nondiagonal components 
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where 



as well as letting a 7^ a. It is not easy to estimate the effect of these terms. For stationary 
configurations, the Einstein equations for the nondiagonal metric components are written 
as follows (see e.g. p4| ): 

^- = 16vr/.v/^^-^. , (52) 

a = det (7ij , (53) 

(Tij = Qij + hgiQj , (54) 

^i = — - , (55) 

9tt 

h = -gtt , (56) 

w'' = a'^a^Wi , (57) 

= ^ _ ^ (58) 

It* 

V^(l — 5'jM-') 
i;2 = aijvV . (60) 

These equations can be considered to be linear equations for gi with source terms propor- 
tional to the velocity f * if h and ctjj are assumed to be known. This implies that even when 
matter velocity is along the azimuthal direction, i.e. f"^ ^ 0,v^ = v^ = 0, other nondiagonal 
components of the metric will emerge. Thus, in general, three nondiagonal components may 
be of the same order because they are proportional to the value of the 3-velocity. However, 
since effect of the nondiagonal components is physically interpreted as dragging of the iner- 
tial frame to the corresponding direction, it may be natural to consider that the nondiagonal 
component along the direction of the matter velocity dominates over nondiagonal compo- 
nents along two other directions. Quantitative estimation should be done by extending the 
present scheme or devising new schemes. 

Concerning the form of the metric, we need to discuss the relation between our metric 
and the CFC scheme. In the investigations which employ the CFC, the following form of 
the metric has been used pll , |33| : 
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ds'^ = -a^dt^ + -iij{dx' - Jdt) {dx^ - Jdt) , (61) 

where 

%, = ^V., , (62) 

and Latin indices run from 1 to 3. Here fij is the flat space metric and \1/ is the conformal 
factor, which is independent of t. The spatial part of this metric is written as: 

g.jdx'dx^ = ^^(cir^ + r^de^ + r^ sin^ ed'f^) . (63) 

On the other hand, the spatial part of our metric is: 

gijdx'dx^ = e^^dr^ + e^'^r^dO^ + e'^^r^ sin^ Odip^ . (64) 

Therefore by comparing a and (3 obtained from our code, we can see how the CFG is satisfied. 
Cook et al. [^ did the same estimate for axisymmetric models and found that the CFC is 
well satisfied for rapidly rotating and general relativistic stars. 

In Figure |^ distributions of the metric functions in the meridional plane for a selected 
model are shown. In Figure ^ relative difference of | exp(2a) — exp(2/?)|/exp(2a) is plotted 
against the distance from the rotation axis. As seen from this figure, distributions of two 
metric components a and (3 are not similar. As far as k < 0.3, relative differences between 
a and (3 are within 5% not only for axisymmetric stars but also for binary star systems. 
Thus in this range of gravity strength, the CFC may be a good approximation. However, 
for models with n ~ 0.5, deviations are larger than 10%. We note that this conclusion is 
obtained for our metric (^ which is not the most general one. However, our results show 
that one should be careful when the CFC is used for extremely relativistic models. 

The local flatness near the rotational axis requires that 

a = P , on the rotation axis. (65) 

From Figures |^ and ^ the value of a coincides with that of /3 to within less than 2% on 
the rotational axis. Thus our numerical code satisfies the local flatness condition to this 
accuracy. 
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As seen from Equations (p3D-(Pq), we must cover the whole space from the origin to 
infinity in the integrations. In actual computations, however, we have covered only a finite 
region, i.e. the region of integration is limited to < f < 2. To check the error due to this 
truncation of the region, we have solved a few models by setting the size of the region to be 
twice as large as the former cases : < f < 4. The difference of nondimensional physical 
quantities, such as the mass, the angular momentum and the angular velocity, between 
these two cases are less than 0.5%. Concerning the approximate angular momentum, since 
its value depends on the distributions of the metric functions as seen from Equation (^8|), 
the obtained approximate angular momentum contains some other error. When we change 



the integral region from 0<f<2to0<f<lin (pSf ), difference is less than 0.2%. Our 
obtained values can be considered to be accurate enough. 

As for the equation of state, we have chosen the relation (Q) as our polytropes. Quanti- 
tative values in this paper are different from those obtained from other choices of polytrope 
such as: 

e = p, + Np, p = K'pI-"^"" . (66) 

As shown in the previous section, although there exists a turning point along the 
quasiequilibrium sequence for stiff equations of state, no turning points appear for softer 
equations of state. The existence of turning points along properly chosen equilibrium se- 
quences are deeply related to the occurrence of some kind of instability (see e.g. [0,Q). 
Therefore, it is important to develop numerical schemes which can compute unstable equi- 
librium configurations. In this sense, our present scheme is the one which is extended to 
handle realistic equilibrium sequences. 

B. Conclusion 

In this paper we have presented a new numerical scheme to handle 3D configurations in 
quasiequilibrium states in general relativity. By using the new scheme, we have succeeded 



24 



in obtaining quasiequilibrium sequences of synchronously rotating binary star systems. We 
have treated restricted situations about 1) the metric, 2) the equation of state, and 3) the 
velocity field for the binary systems. The next step of our investigations is to remove some or 
all of these restrictions by extending the formulation to the more general form of the metric, 
including realistic equations of state and/or treating irrotational binary star systems. At 
that stage, we will be able to study turning points or critical configurations by applying the 
extended formulation. 
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APPENDIX A: RICCI TENSOR AND EINSTEIN EQUATIONS IN THE 

TETRAD SYSTEM 



Ricci tensors associated with the metric (0^ 
system as follows: 

2dv 



-R(C)(»y)5 is explicitly written in the tetrad 
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where we have used the orthonormal tetrad defined as: 

' e-^ 0, 0, (cu-fi)e-'^ 

0, e"", 0, 

0, 0, ^, 

0, 0, 0, ^ 
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Any tensor, /i^^,, is transformed to that in the tetrad system as follows: 



^(a)(/3) - ht;V^{a)^{f3) 



(A13) 
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The energy-momentum tensor for the quasiequihbrium configurations in this paper are: 

^ + ^"' (AM) 

(A15) 
(A16) 

(A17) 
(A18) 

Einstein equation is formally written as: 

R{Oiv) = 8^(^(C)W - 2^iOiv)T) , (A19) 

where T = T.5n . We have used the following combinations to get elliptic type equations 





^ww - 


l-i;2 ' 


T(r)(r) 


= T(e)[e) = 


= P , 




^{'p)('p) = 




T(t){f) 


= ^{^)W = 


{e + p)v 




others = 


= . 



(0)-(|g) and (0). 

Rim + Rm^) = 87r(£ + P) Y3^ . (A20) 

^(t)(t) - R^M = 167rp , (A21) 

R{t)(t) + -R(r)(r) + R(e){e) — R(^){ip) = 8TT{e + p) , (A22) 

i?(,)(^)=87r(£ + p)^^, (A23) 

APPENDIX B: DETAILED NUMERICAL RESULTS FOR BINARY SEQUENCES 



Numerical results of binary sequences are shown in Tables [I- Xlll . 

In these Tables, f^, J, M, Mq, Cl"^ and fx are the dimensionless distance to point A from 
the origin, the dimensionless total angular momentum, the dimensionless total gravitational 
mass, the dimensionless total rest mass, the dimensionless angular velocity and the scale 
parameter, respectively. Here, masses are normalized as: 

M 

^0 = o,, ^\ wo • (B2) 
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FIGURES 
FIG. 1. Schematic view of the binary system. The x— axis is set along the line joining the two 

centers of mass of the two stars and the origin is set on the middle point. The z— axis is along the 

rotation axis. Point A, B are set at the intersections of the surface of the star and the x— axis. The 

inner intersection is point A and the outer one is point B. 

FIG. 2. Dimensionless squared angular velocity ^"^ /{AttGec/c?) is plotted against dimensionless 
angular momentum J /{AttGc^'^Mq /ej )^'^ for the sequences of the axisymmetric stars from 
fp/fe = 0.9375 (nearly spherical) to rp/Vf, = 0.5625 with polytropic indices (a) A^ = 0.0 and 
(b) A^ = 0.5, where rp and r^ are the polar radius and the equatorial radius of axisymmetric 
configurations. The solid line (k = 0.0001), the dotted line {k = 0.25) and the dashed line 
[k = 0.4) are results obtained by using our 3D code. Crosses {k = 0.0001), open circles (k = 0.25) 
and filled circles {k = 0.4) are taken from Komatsu et al. p8| . 

FIG. 3. Same as Figure |2| but for the sequences of Newtonian binary systems from q = 0.0 to 
q = 0.5. The solid line (A^ = 0.0), the dotted line (A^ = 0.5), the short-dashed line (A^ = 1.0) 
and the long-dashed line (A^ = 1.5) are our numerical results. Crosses (A^ = 0.0), filled circles 
(A^ = 0.5), open circles (A^ = 1.0) and filled rectangles (A^ = 1.5) are taken from Hachisu [13|. 



FIG. 4. Same as Figure ^ but for the sequences of relativistic binary systems from q = 0.0 
to q = 0.5 with polytropic indices (a) A' = 0.0, (b) A^ = 0.5 and (c) A^ = 1.0. On each panel, 
K = 0.0001 (Newtonian), 0.1, 0.3 and 0.5 are shown by solid lines. Crosses denote the Newtonian 
results taken from Hachisu |13]. 



FIG. 5. Contours of the energy density on the x — z plane. The units of the distance is r^. (a) 
The contact phase {q = 0.0) and (b) the distant phase {q = 0.5) are shown for A^ = 0.5 polytrope. 
On each panel, contours for k = 0.0001 (Newtonian) and k = 0.5 (relativistic) are shown. 
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FIG. 6. Same as Figure but for A'^ = 0.0 polytropes with different k and different topology, i.e. 
axisymmetric sequences, Jacobi-Dumb-bell shaped sequence and binary sequences. SoHd, dotted 
and dashed curves denote sequences with k = 0.0001 (Newtonian), 0.25 and 0.4, respectively. The 
dash-dotted curve denotes the Jacobi-Dumb-bell shaped sequences taken from Hachisu in 



FIG. 7. (a) Distribution of the metric function exp(2i/) on the equatorial plane is plotted 
against the distance from the rotational axis. Model parameters are A^ = 0.5, n = 0.5 and q = 0.0. 
Different curves correspond to distributions on 99 = 0, vr/12, 27r/12, 37r/12, 47r/12, 57r/12, 67r/12. (b) 
Distribution of the metric function exp(2/3). (c) Distribution of the metric function exp(2a). (d) 
Distribution of the metric function uj/Vt. 

FIG. 8. Distribution of | exp(2a) — exp(2/3)|/exp(2a) on the equatorial plane is plotted against 
the distance from the rotational axis. Model parameters are A^ = 0.5, k = 0.5 and q = 0.0. 
Different curves correspond to distributions on (/? = 0, 37r/12,67r/12. 
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TABLES 
TABLE I. N = 0.0, k = 0.0001 (Newtonian) 



fA J M Mo ^' fx 

l.OOOE-08 1.126E-01 3.246E-06 3.246E-06 2.647E-02 3.984E-04 

3.125E-02 1.121E-01 3.233E-06 3.234E-06 2.644E-02 3.960E-04 

6.250E-02 1.109E-01 3.187E-06 3.187E-06 2.586E-02 3.894E-04 

9.375E-02 1.097E-01 3.137E-06 3.138E-06 2.465E-02 3.836E-04 

1.250E-01 1.082E-01 3.078E-06 3.078E-06 2.275E-02 3.797E-04 

1.562E-01 1.078E-01 3.032E-06 3.033E-06 2.077E-02 3.830E-04 

1.875E-01 1.074E-01 2.991E-06 2.992E-06 1.849E-02 3.893E-04 

2.188E-01 1.078E-01 2.965E-06 2.965E-06 1.631E-02 4.013E-04 

2.500E-01 1.087E-01 2.938E-06 2.938E-06 1.423E-02 4.176E-04 

2.812E-01 1.096E-01 2.905E-06 2.906E-06 1.213E-02 4.379E-04 

3.125E-01 1.112E-01 2.874E-06 2.874E-06 1.033E-02 4.620E-04 

3.438E-01 1.140E-01 2.825E-06 2.826E-06 8.790E-03 4.914E-04 

3.750E-01 1.150E-01 2.780E-06 2.780E-06 7.110E-03 5.261E-04 

4.062E-01 1.182E-01 2.752E-06 2.753E-06 5.880E-03 5.713E-04 

4.375E-01 1.212E-01 2.750E-06 2.750E-06 4.808E-03 6.262E-04 

4.688E-01 1.201E-01 2.780E-06 2.781E-06 3.535E-03 7.014E-04 

5.000E-01 1.295E-01 2.820E-06 2.821E-06 3.089E-03 7.916E-04 
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TABLE II. N = 0.5, k = 0.0001 (Newtonian) 



TA J M Mo ^2 f^ 

l.OOOE-08 1.145E-01 4.342E-06 4.342E-06 1.668E-02 6.474E-04 

3.125E-02 1.142E-01 4.330E-06 4.330E-06 1.647E-02 6.464E-04 

6.250E-02 1.137E-01 4.297E-06 4.298E-06 1.582E-02 6.455E-04 

9.375E-02 1.132E-01 4.256E-06 4.257E-06 1.480E-02 6.477E-04 

1.250E-01 1.130E-01 4.209E-06 4.210E-06 1.353E-02 6.549E-04 

1.562E-01 1.131E-01 4.168E-06 4.169E-06 1.213E-02 6.684E-04 

1.875E-01 1.136E-01 4.133E-06 4.134E-06 1.070E-02 6.891E-04 

2.188E-01 1.145E-01 4.101E-06 4.102E-06 9.309E-03 7.168E-04 

2.500E-01 1.158E-01 4.071E-06 4.072E-06 7.998E-03 7.523E-04 

2.812E-01 1.175E-01 4.041E-06 4.042E-06 6.800E-03 7.958E-04 

3.125E-01 1.195E-01 4.009E-06 4.010E-06 5.711E-03 8.484E-04 

3.438E-01 1.219E-01 3.980E-06 3.980E-06 4.746E-03 9.116E-04 

3.750E-01 1.246E-01 3.957E-06 3.957E-06 3.900E-03 9.877E-04 

4.062E-01 1.278E-01 3.948E-06 3.949E-06 3.188E-03 1.079E-03 

4.375E-01 1.316E-01 3.955E-06 3.955E-06 2.592E-03 1.190E-03 

4.688E-01 1.349E-01 3.968E-06 3.969E-06 2.052E-03 1.325E-03 

5.000E-01 1.408E-01 3.982E-06 3.982E-06 1.666E-03 1.486E-03 
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TABLE III. N = 1.0, K = 0.0001 (Newtonian) 



TA J M Mo ^2 f^ 

l.OOOE-08 1.182E-01 5.336E-06 5.337E-06 9.957E-03 1.028E-03 

3.125E-02 1.182E-01 5.327E-06 5.328E-06 9.802E-03 1.030E-03 

6.250E-02 1.183E-01 5.303E-06 5.304E-06 9.359E-03 1.037E-03 

9.375E-02 1.184E-01 5.274E-06 5.274E-06 8.702E-03 1.051E-03 

1.250E-01 1.189E-01 5.238E-06 5.239E-06 7.903E-03 1.073E-03 

1.562E-01 1.197E-01 5.206E-06 5.207E-06 7.044E-03 1.105E-03 

1.875E-01 1.208E-01 5.178E-06 5.178E-06 6.181E-03 1.147E-03 

2.188E-01 1.223E-01 5.148E-06 5.149E-06 5.345E-03 1.201E-03 

2.500E-01 1.241E-01 5.121E-06 5.122E-06 4.568E-03 1.268E-03 

2.812E-01 1.264E-01 5.096E-06 5.097E-06 3.861E-03 1.349E-03 

3.125E-01 1.289E-01 5.073E-06 5.074E-06 3.231E-03 1.446E-03 

3.438E-01 1.319E-01 5.055E-06 5.056E-06 2.679E-03 1.561E-03 

3.750E-01 1.352E-01 5.043E-06 5.043E-06 2.201E-03 1.697E-03 

4.062E-01 1.390E-01 5.037E-06 5.038E-06 1.796E-03 1.859E-03 

4.375E-01 1.431E-01 5.035E-06 5.036E-06 1.451E-03 2.052E-03 

4.688E-01 1.477E-01 5.031E-06 5.032E-06 1.158E-03 2.281E-03 

5.000E-01 1.532E-01 5.021E-06 5.022E-06 9.183E-04 2.558E-03 
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TABLE IV. N = 1.5, n = 0.0001 (Newtonian) 



TA J M Mo ^2 f^ 

l.OOOE-08 1.246E-01 6.282E-06 6.282E-06 5.596E-03 1.662E-03 

3.125E-02 1.247E-01 6.278E-06 6.278E-06 5.509E-03 1.667E-03 

6.250E-02 1.251E-01 6.260E-06 6.261E-06 5.252E-03 1.684E-03 

9.375E-02 1.258E-01 6.240E-06 6.240E-06 4.878E-03 1.714E-03 

1.250E-01 1.268E-01 6.214E-06 6.215E-06 4.422E-03 1.758E-03 

1.562E-01 1.281E-01 6.191E-06 6.192E-06 3.936E-03 1.818E-03 

1.875E-01 1.297E-01 6.170E-06 6.170E-06 3.448E-03 1.896E-03 

2.188E-01 1.318E-01 6.146E-06 6.146E-06 2.975E-03 1.993E-03 

2.500E-01 1.342E-01 6.124E-06 6.125E-06 2.536E-03 2.110E-03 

2.812E-01 1.370E-01 6.105E-06 6.106E-06 2.140E-03 2.252E-03 

3.125E-01 1.401E-01 6.089E-06 6.089E-06 1.788E-03 2.419E-03 

3.438E-01 1.437E-01 6.075E-06 6.076E-06 1.481E-03 2.617E-03 

3.750E-01 1.476E-01 6.064E-06 6.065E-06 1.215E-03 2.850E-03 

4.062E-01 1.519E-01 6.055E-06 6.056E-06 9.891E-04 3.125E-03 

4.375E-01 1.567E-01 6.046E-06 6.046E-06 7.965E-04 3.453E-03 

4.688E-01 1.621E-01 6.036E-06 6.036E-06 6.353E-04 3.844E-03 

5.000E-01 1.680E-01 6.024E-06 6.024E-06 5.006E-04 4.314E-03 
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TABLE V. iV = 0.0, k = 0.1 



VA J M Mo ^' fx 

l.OOOE-08 1.316E-01 5.411E-02 5.997E-02 2.794E-02 1.531E-01 

3.125E-02 1.307E-01 5.386E-02 5.968E-02 2.742E-02 1.533E-01 

6.250E-02 1.296E-01 5.371E-02 5.953E-02 2.716E-02 1.519E-01 

9.375E-02 1.286E-01 5.317E-02 5.892E-02 2.553E-02 1.528E-01 

1.250E-01 1.274E-01 5.263E-02 5.829E-02 2.375E-02 1.543E-01 

1.562E-01 1.268E-01 5.215E-02 5.774E-02 2.145E-02 1.579E-01 

1.875E-01 1.270E-01 5.177E-02 5.728E-02 1.912E-02 1.634E-01 

2.188E-01 1.280E-01 5.153E-02 5.697E-02 1.687E-02 1.712E-01 

2.500E-01 1.291E-01 5.132E-02 5.672E-02 1.463E-02 1.804E-01 

2.812E-01 1.306E-01 5.117E-02 5.653E-02 1.250E-02 1.922E-01 

3.125E-01 1.325E-01 5.089E-02 5.619E-02 1.056E-02 2.060E-01 

3.438E-01 1.351E-01 5.056E-02 5.578E-02 8.922E-03 2.218E-01 

3.750E-01 1.375E-01 5.012E-02 5.526E-02 7.359E-03 2.415E-01 

4.062E-01 1.405E-01 4.981E-02 5.489E-02 5.996E-03 2.653E-01 

4.375E-01 1.439E-01 4.991E-02 5.500E-02 4.883E-03 2.929E-01 

4.688E-01 1.435E-01 5.020E-02 5.543E-02 3.631E-03 3.290E-01 

5.000E-01 1.556E-01 5.070E-02 5.584E-02 3.180E-03 3.713E-01 
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TABLE VI. iV = 0.5, k = 0.1 



VA J M Mo ^' fx 

l.OOOE-08 1.438E-01 7.006E-02 7.664E-02 1.631E-02 2.457E-01 

3.125E-02 1.436E-01 6.998E-02 7.655E-02 1.608E-02 2.464E-01 

6.250E-02 1.432E-01 6.975E-02 7.629E-02 1.547E-02 2.486E-01 

9.375E-02 1.431E-01 6.943E-02 7.593E-02 1.443E-02 2.533E-01 

1.250E-01 1.432E-01 6.912E-02 7.557E-02 1.318E-02 2.602E-01 

1.562E-01 1.437E-01 6.882E-02 7.523E-02 1.180E-02 2.699E-01 

1.875E-01 1.448E-01 6.859E-02 7.495E-02 1.040E-02 2.826E-01 

2.188E-01 1.462E-01 6.842E-02 7.475E-02 9.029E-03 2.983E-01 

2.500E-01 1.481E-01 6.831E-02 7.461E-02 7.750E-03 3.174E-01 

2.812E-01 1.505E-01 6.819E-02 7.446E-02 6.574E-03 3.404E-01 

3.125E-01 1.533E-01 6.805E-02 7.427E-02 5.518E-03 3.679E-01 

3.438E-01 1.564E-01 6.791E-02 7.411E-02 4.582E-03 4.003E-01 

3.750E-01 1.598E-01 6.783E-02 7.399E-02 3.760E-03 4.386E-01 

4.062E-01 1.640E-01 6.784E-02 7.399E-02 3.068E-03 4.840E-01 

4.375E-01 1.689E-01 6.799E-02 7.414E-02 2.487E-03 5.376E-01 

4.688E-01 1.731E-01 6.820E-02 7.439E-02 1.962E-03 6.022E-01 

5.000E-01 1.807E-01 6.842E-02 7.459E-02 1.588E-03 6.787E-01 
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TABLE VII. N = 1.0, k = 0.1 



fA J M Mo ^2 f^ 

l.OOOE-08 1.612E-01 8.428E-02 9.030E-02 8.716E-03 4.072E-01 

3.125E-02 1.631E-01 8.461E-02 9.060E-02 8.861E-03 4.080E-01 

6.250E-02 1.618E-01 8.435E-02 9.037E-02 8.276E-03 4.151E-01 

9.375E-02 1.624E-01 8.421E-02 9.022E-02 7.712E-03 4.252E-01 

1.250E-01 1.634E-01 8.407E-02 9.005E-02 7.022E-03 4.397E-01 

1.562E-01 1.649E-01 8.395E-02 8.992E-02 6.269E-03 4.588E-01 

1.875E-01 1.668E-01 8.384E-02 8.979E-02 5.503E-03 4.829E-01 

2.188E-01 1.691E-01 8.378E-02 8.971E-02 4.763E-03 5.122E-01 

2.500E-01 1.719E-01 8.375E-02 8.967E-02 4.072E-03 5.474E-01 

2.812E-01 1.751E-01 8.375E-02 8.966E-02 3.443E-03 5.891E-01 

3.125E-01 1.788E-01 8.375E-02 8.965E-02 2.879E-03 6.386E-01 

3.438E-01 1.830E-01 8.376E-02 8.965E-02 2.385E-03 6.969E-01 

3.750E-01 1.876E-01 8.381E-02 8.969E-02 1.957E-03 7.653E-01 

4.062E-01 1.929E-01 8.391E-02 8.979E-02 1.593E-03 8.456E-01 

4.375E-01 1.986E-01 8.404E-02 8.992E-02 1.284E-03 9.405E-01 

4.688E-01 2.048E-01 8.417E-02 9.007E-02 1.023E-03 1.054E+00 

5.000E-01 2.119E-01 8.429E-02 9.018E-02 8.073E-04 1.190E+00 



41 



TABLE VIII. N = 0.0, k = 0.3 



fA J M Mo ^' fx 

l.OOOE-08 1.638E-01 1.231E-01 1.469E-01 2.870E-02 1.447E-01 

3.125E-02 1.629E-01 1.227E-01 1.466E-01 2.825E-02 1.452E-01 

6.250E-02 1.620E-01 1.223E-01 1.461E-01 2.773E-02 1.457E-01 

9.375E-02 1.613E-01 1.217E-01 1.455E-01 2.604E-02 1.491E-01 

1.250E-01 1.600E-01 1.212E-01 1.450E-01 2.406E-02 1.536E-01 

1.562E-01 1.598E-01 1.209E-01 1.448E-01 2.170E-02 1.600E-01 

1.875E-01 1.603E-01 1.206E-01 1.445E-01 1.919E-02 1.688E-01 

2.188E-01 1.609E-01 1.206E-01 1.448E-01 1.674E-02 1.790E-01 

2.500E-01 1.627E-01 1.208E-01 1.450E-01 1.441E-02 1.924E-01 

2.812E-01 1.651E-01 1.210E-01 1.453E-01 1.240E-02 2.073E-01 

3.125E-01 1.673E-01 1.210E-01 1.455E-01 1.040E-02 2.260E-01 

3.438E-01 1.697E-01 1.208E-01 1.453E-01 8.631E-03 2.479E-01 

3.750E-01 1.725E-01 1.206E-01 1.451E-01 7.139E-03 2.728E-01 

4.062E-01 1.755E-01 1.204E-01 1.450E-01 5.749E-03 3.048E-01 

4.375E-01 1.778E-01 1.205E-01 1.454E-01 4.562E-03 3.405E-01 

4.688E-01 1.801E-01 1.208E-01 1.462E-01 3.489E-03 3.861E-01 

5.000E-01 1.907E-01 1.220E-01 1.474E-01 2.910E-03 4.366E-01 
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TABLE IX. N = 0.5, k = 0.3 



fA J M Mo ^' fx 

l.OOOE-08 2.013E-01 1.556E-01 1.762E-01 1.429E-02 2.520E-01 

3.125E-02 2.013E-01 1.555E-01 1.762E-01 1.410E-02 2.536E-01 

6.250E-02 2.013E-01 1.554E-01 1.761E-01 1.359E-02 2.577E-01 

9.375E-02 2.015E-01 1.554E-01 1.762E-01 1.269E-02 2.656E-01 

1.250E-01 2.021E-01 1.554E-01 1.765E-01 1.159E-02 2.765E-01 

1.562E-01 2.031E-01 1.555E-01 1.768E-01 1.037E-02 2.911E-01 

1.875E-01 2.048E-01 1.556E-01 1.773E-01 9.128E-03 3.090E-01 

2.188E-01 2.069E-01 1.559E-01 1.778E-01 7.901E-03 3.308E-01 

2.500E-01 2.096E-01 1.562E-01 1.784E-01 6.760E-03 3.569E-01 

2.812E-01 2.127E-01 1.565E-01 1.790E-01 5.705E-03 3.878E-01 

3.125E-01 2.163E-01 1.568E-01 1.796E-01 4.770E-03 4.243E-01 

3.438E-01 2.203E-01 1.571E-01 1.801E-01 3.942E-03 4.673E-01 

3.750E-01 2.247E-01 1.574E-01 1.808E-01 3.220E-03 5.179E-01 

4.062E-01 2.298E-01 1.577E-01 1.814E-01 2.607E-03 5.773E-01 

4.375E-01 2.359E-01 1.582E-01 1.822E-01 2.096E-03 6.472E-01 

4.688E-01 2.413E-01 1.587E-01 1.830E-01 1.646E-03 7.308E-01 

5.000E-01 2.499E-01 1.594E-01 1.839E-01 1.308E-03 8.297E-01 
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TABLE X. iV = 1.0, K = 0.3 



VA J M Mo ^2 f^ 

l.OOOE-08 2.557E-01 1.882E-01 2.007E-01 6.089E-03 4.906E-01 

3.125E-02 2.594E-01 1.881E-01 1.999E-01 6.206E-03 4.901E-01 

6.250E-02 2.576E-01 1.882E-01 2.007E-01 5.827E-03 5.021E-01 

9.375E-02 2.589E-01 1.884E-01 2.011E-01 5.453E-03 5.174E-01 

1.250E-01 2.606E-01 1.886E-01 2.016E-01 4.984E-03 5.395E-01 

1.562E-01 2.630E-01 1.889E-01 2.022E-01 4.466E-03 5.681E-01 

1.875E-01 2.660E-01 1.893E-01 2.029E-01 3.930E-03 6.036E-01 

2.188E-01 2.697E-01 1.897E-01 2.036E-01 3.407E-03 6.457E-01 

2.500E-01 2.738E-01 1.902E-01 2.045E-01 2.908E-03 6.975E-01 

2.812E-01 2.784E-01 1.906E-01 2.053E-01 2.452E-03 7.583E-01 

3.125E-01 2.837E-01 1.911E-01 2.062E-01 2.044E-03 8.292E-01 

3.438E-01 2.897E-01 1.916E-01 2.071E-01 1.687E-03 9.117E-01 

3.750E-01 2.966E-01 1.922E-01 2.080E-01 1.379E-03 1.008E+00 

4.062E-01 3.041E-01 1.928E-01 2.089E-01 1.117E-03 1.120E+00 

4.375E-01 3.125E-01 1.934E-01 2.099E-01 8.958E-04 1.252E+00 

4.688E-01 3.214E-01 1.940E-01 2.109E-01 7.092E-04 1.408E+00 

5.000E-01 3.314E-01 1.947E-01 2.119E-01 5.562E-04 1.590E+00 
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TABLE XL N = 0.0, n = 0.5 (Highly relativistic) 



rA 



J 



M 



Mo 



n^ 



Jx 



LOOOE-08 
3.125E-02 
6.250E-02 
9.375E-02 
1.250E-01 
1.562E-01 
L875E-01 
2.188E-01 
2.500E-01 
2.812E-01 
3.125E-01 
3.438E-01 
3.750E-01 
4.062E-01 
4.375E-01 
4.688E-01 
5.000E-01 



1.911E-01 
1.901E-01 
1.889E-01 
1.880E-01 
1.871E-01 
1.867E-01 
1.871E-01 
1.883E-01 
1.900E-01 
1.922E-01 
1.952E-01 
1.976E-01 
2.007E-01 
2.029E-01 
2.060E-01 
2.081E-01 
2.187E-01 



1.523E-01 
1.519E-01 
1.516E-01 
1.516E-01 
1.514E-01 
1.516E-01 
1.518E-01 
1.523E-01 
1.528E-01 
1.534E-01 
L538E-01 
L539E-01 
L539E-01 
L537E-01 
1.541E-01 
1.548E-01 
L559E-01 



1.833E-01 
1.831E-01 
1.829E-01 
L836E-01 
L837E-01 
1.846E-01 
1.854E-01 
1.865E-01 
1.877E-01 
1.890E-01 
1.898E-01 
1.903E-01 
1.908E-01 
1.912E-01 
1.922E-01 
1.940E-01 
1.951E-01 



2.866E-02 
2.819E-02 
2.749E-02 
2.564E-02 
2.362E-02 
2.122E-02 
L875E-02 
1.623E-02 
1.396E-02 
1.183E-02 
9.972E-03 
8.243E-03 
6.785E-03 
5.397E-03 
4.280E-03 
3.310E-03 
2.703E-03 



1.209E-01 
1.214E-01 
1.225E-01 
1.263E-01 
1.315E-01 
1.385E-01 
1.473E-01 
1.586E-01 
1.709E-01 
1.866E-01 
2.049E-01 
2.270E-01 
2.520E-01 
2.835E-01 
3.191E-01 
3.618E-01 
4.115E-01 
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TABLE XII. N = 0.5, n = 0.5 (Highly relativistic) 



rA 



J 



M 



Mo 



n^ 



Jx 



l.OOOE-08 
3.125E-02 
6.250E-02 
9.375E-02 
1.250E-01 
1.562E-01 
1.875E-01 
2.188E-01 
2.500E-01 
2.812E-01 
3.125E-01 
3.438E-01 
3.750E-01 
4.062E-01 
4.375E-01 
4.688E-01 
5.000E-01 



2.589E-01 
2.586E-01 
2.586E-01 
2.590E-01 
2.598E-01 
2.613E-01 
2.633E-01 
2.660E-01 
2.694E-01 
2.728E-01 
2.769E-01 
2.814E-01 
2.861E-01 
2.920E-01 
2.990E-01 
3.051E-01 
3.141E-01 



1.956E-01 
1.955E-01 
1.956E-01 
1.960E-01 
1.965E-01 
1.970E-01 
1.977E-01 
1.984E-01 
1.992E-01 
2.000E-01 
2.007E-01 
2.014E-01 
2.022E-01 
2.029E-01 
2.038E-01 
2.045E-01 
2.055E-01 



2.161E-01 
2.161E-01 
2.166E-01 
2.175E-01 
2.187E-01 
2.200E-01 
2.215E-01 
2.232E-01 
2.249E-01 
2.266E-01 
2.283E-01 
2.299E-01 
2.317E-01 
2.333E-01 
2.350E-01 
2.367E-01 
2.385E-01 



1.217E-02 
1.199E-02 
1.154E-02 
1.079E-02 
9.869E-03 
8.839E-03 
7.779E-03 
6.733E-03 
5.750E-03 
4.839E-03 
4.029E-03 
3.313E-03 
2.688E-03 
2.164E-03 
1.730E-03 
1.350E-03 
1.058E-03 



2.379E-01 
2.394E-01 
2.439E-01 
2.524E-01 
2.640E-01 
2.791E-01 
2.981E-01 
3.209E-01 
3.480E-01 
3.805E-01 
4.189E-01 
4.641E-01 
5.173E-01 
5.790E-01 
6.517E-01 
7.385E-01 
8.408E-01 
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TABLE XIII. N = 1.0, k = 0.5 (Highly relativistic) 



rA 



J 



M 



Mo 



n^ 



Jx 



l.OOOE-08 
3.125E-02 
6.250E-02 
9.375E-02 
1.250E-01 
1.562E-01 
1.875E-01 
2.188E-01 
2.500E-01 
2.812E-01 
3.125E-01 
3.438E-01 
3.750E-01 
4.062E-01 
4.375E-01 
4.688E-01 
5.000E-01 



3.578E-01 
3.583E-01 
3.592E-01 
3.608E-01 
3.637E-01 
3.667E-01 
3.711E-01 
3.757E-01 
3.816E-01 
3.878E-01 
3.927E-01 
4.025E-01 
4.120E-01 
4.232E-01 
4.357E-01 
4.489E-01 
4.652E-01 



2.399E-01 
2.403E-01 
2.406E-01 
2.410E-01 
2.415E-01 
2.422E-01 
2.430E-01 
2.438E-01 
2.446E-01 
2.455E-01 
2.466E-01 
2.473E-01 
2.481E-01 
2.489E-01 
2.498E-01 
2.507E-01 
2.517E-01 



2.437E-01 
2.440E-01 
2.446E-01 
2.455E-01 
2.465E-01 
2.479E-01 
2.494E-01 
2.511E-01 
2.528E-01 
2.546E-01 
2.569E-01 
2.582E-01 
2.598E-01 
2.613E-01 
2.628E-01 
2.644E-01 
2.659E-01 



4.250E-03 
4.206E-03 
4.048E-03 
3.800E-03 
3.489E-03 
3.136E-03 
2.769E-03 
2.402E-03 
2.054E-03 
1.730E-03 
1.424E-03 
1.185E-03 
9.668E-04 
7.820E-04 
6.269E-04 
4.971E-04 
3.926E-04 



5.513E-01 

5.555E-01 

5.664E-01 

5.844E-01 

6.089E-01 

6.422E-01 

6.824E-01 

7.325E-01 

7.906E-01 

8.601E-01 

9.449E-01 

1.037E+00 

1.145E+00 

1.268E+00 

1.410E+00 

1.577E+00 

1.767E+00 
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